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Abstract 

In recent work, we have developed a variational principle for large N 
multi-matrix models based on the extremization of non-commutative en- 
tropy. Here, we test the simplest variational ansatz for our entropic varia- 
tional principle with Monte-Carlo measurements. In particular, we study 
the two matrix model with action tr {A\ + A\) — i [A\ , A2] 2 ] which has 
not been exactly solved. We estimate the expectation values of traces of 
products of matrices and also those of traces of products of exponentials 
of matrices (Wilson loop operators) . These are compared with a Monte- 
Carlo simulation. We find that the simplest wignerian variational ansatz 
provides a remarkably good estimate for observables when m 2 is of order 
unity or more. For small values of m 2 the wignerian ansatz is not a good 
approximation: the measured correlations grow without bound, reflect- 
ing the non-convergence of matrix integrals defining the pure commutator 
squared action. Comparison of this ansatz with the exact solution of a two 
matrix model studied by Mehta is also summarized. Here the wignerian 
ansatz is a good approximation both for strong and weak coupling. 

KEYWORDS: multi matrix models, Yang-Mills integrals, Yang-Mills the- 
ory, M(atrix) theory, large N limit, entropy, variational principle, Monte-Carlo 
integration. 
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1 Introduction 



Matrix models are of long standing interest in several branches of physics and 
mathematics. 

Early work of Wigner, Dyson and others introduced random matrices in the 
study of statistical properties of highly excited energy levels of nuclei pQ . 

In the early 1970s the work of 't Hooft HO] on the large N limit of QCD (N 
is the number of colors, gluon fields are NxN matrices), made the study of large 
N matrix field theories a central theme in understanding the non-perturbative 
dynamics of non-abelian gauge theories. This also gave the first indication that 
it is in the large N limit that a gauge theory may have the long sought dual 
description as a string theory. 

Important progress was made in the late 1970s and early 1980s stemming 
from Migdal and Makeenko's work on the factorized loop equations of large N 
gauge theories 0]. The work of Sakita and Jevicki on the collective field 
formalism of large N field theories and that of Cvitanovic and collaborators 
[S] also bears mention from this period. Eguchi and Kawai's proposal on 
reducing a matrix field theory to a matrix model with a finite number of degrees 
of freedom but in the large N limit has been a recurring theme ever since. 

Important breakthroughs in the study of random surfaces, two dimensional 
string theory and two dimensional gravity coupled to matter were made in the 
late 80s and early 1990s (see [H] for a review) The planar Feynman graph ex- 
pansion of large N matrix models was used as a way of descretizing a two 
dimensional surface. Models with one or a finite number of matrices and the 
c = 1 quantum mechanics of a single matrix were of importance in these develop- 
ments. In addition to the large N limit, the double scaling limit was developed 
to study the surfaces obtained in the continuum limit. 

From the early 1990s onwards, the work of mathematicians including Voiculescu 
and collaborators on von Neumann algebras lead to the development of the field 
of non-commutative probability theory [§]. Large N matrix models are natu- 
ral examples of non-commutative probability theories. Random matrices also 
have deep connections to the statistical properties of zeros of the Riemann zeta 
function pQ. 

In the mid 1990s, supersymmetric matrix models we proposed as non-perturbative 
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definitions of M-theory and superstring theory El E] • Bosonic matrix models 
are also studied El EI El El El E] as a first step towards understanding 
these supersymmetric matrix models. 

Recently, interest in matrix models has been revived in several contexts. 

Work of Dijkgraaf and Vafa |18U19) has shown that the effective superpoten- 
tial of M = 1 supersymmetric gauge theories with an adjoint chiral superficld 
can be calculated exactly from a bosonic one matrix model in the large TV limit. 
This correspondence also has extensions to multi-matrix models, for example 
the glueball superpotential of Af = 1* gauge theory was determined from the 
partition function of the corresponding three-matrix model |19U22ll2Tj . Bosonic 
multi-matrix models also arise in the context of quiver gauge theories |2()| . 

The c — 1 matrix model has been revived in the context of two dimensional 
super-gravity coupled to c = 1 matter |23|. c = 1 matrix models are also of 
current interest in studying unstable D-branes and the phenomenon of tachyon 
condensation |24l I25j . 

We may also regard finite matrix models as zero momentum limits of matrix 
field theories such as Yang-Mills theory. They provide a testing ground for new 
ideas on matrix field theories. 

While matrix models have been extending their influence into a wide variety 
of contexts, there have been attempts to solve matrix models exactly. For 
instance, the work of Brezin et.al. j2EJ, Mehta |2Z1, Kazakov Staudacher 
|28j and others has shed light on the partition functions and certain special 
classes of correlations in one and two matrix models. The two matrix model 
studied by Mehta is actually a special case of a class of chain multi-matrix 
models which can all be solved exactly, leading in the continuum limit to the 
solution of the c = 1 matrix quantum mechanics found by Brezin et. al. 

However, it has proven difficult to determine the correlations of a generic 
multi-matrix model analytically. In the 't Hooft limit the fluctuations in invari- 
ant observables becomes small and the theory becomes classical while retain- 
ing the quantum fluctuations in h and certain other non-perturbative features. 
Aside from some special cases where exact solutions for the partition function 
and certain correlations have been possible, multi-matrix models, even in the 
large N limit are not exactly solved. This should not be surprising, since they 
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represent complicated classical dynamical systems. Monte-carlo integration is 
probably the best available means of obtaining numerically accurate predictions. 
Their drawback is that some of the mathematical structures are not revealed. It 
would be useful to have a middle ground: an approximation method that pro- 
vides both qualitative and mathematical insights and a quantitative estimate 
for correlations. 

In a previous paper (123, see also we derived a variational prin- 

ciple for large N multi-matrix models. This involved many new theoretical 
ideas, especially the role played by the automorphism group of the free algebra 
generated by the matrices. The most important physical principle learnt was 
that the cohomology of this group is a non-commutative entropy. It was the 
crucial element in deriving a classical action for large N matrix models. This 
allowed us to obtain variational approximations for the correlation tensors of 
multi- matrix models in the large TV limit. Thus we have a self contained formu- 
lation and method of approximate solution for matrix models in the N — » oo 
limit. It deals directly with correlations, rather than matrix elements. There is 
no more any reference to integrations over matrices nor the principle of unitary 
invariance. Rather, we have a classical theory on the configuration space of 
non-commutative probability distributions. The correlation tensors are coordi- 
nates on this space. They are determined by a constrained maximization of a 
non-commutative analogue of entropy. The action of the original matrix model 
is encoded in the constraints. 

To solve this extremization problem approximately, we maximize the en- 
tropy on a conveniently chosen finite-dimensional sub-space of the configuration 
space. The simplest such choice is the subspace corresponding to the wignerian 
correlations. 

In |30j we compared our variational ansatz with the exact solution of a two 
matrix model studied by Mehta |27| . The variational ansatz gave a reasonably 
good approximation both for strong and weak coupling (summarized in section 
[5J. This gives us the confidence to test our ansatz for models that have not 
been exactly solved. In this paper we quantitatively test the wignerian vari- 
ational ansatz with correlations measured using Monte-Carlo integration for a 
bosonic two matrix model with action tr \^m 2 (A\ + A?,) — j\A\,A2\ 2 }. We 
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study a two matrix model since it is the simplest multi-matrix model. We pick 
this action since it mimics that of the zero momentum limit of Yang-Mills the- 
ory. Moreover, this is the simplest two matrix model that is not exactly solved 
and also shares the derivation property of Yang-Mills theory [331 • We use 
the Metropolis algorithm implemented on a personal computer to test our new 
variational principle with the simplest of ansatze. We compare measured and 
variational two and four point correlations (eg. < -S^^li^^b A\ >) and also the 
expectation values of Wilson loop operators (eg. < Jig*' *e 2 e — *e - 2 >) 
in the large N limit. We find that it works remarkably well for m 2 of order 1 or 
more. As m 2 — > 0, the measured correlations appear to grow without bound, re- 
flecting the divergence in the matrix integrals that define the pure commutator 
squared model. This is not captured by the wignerian ansatz. 

We now summarize the framework we use to study large N matrix models. 



2 Large N Matrix Models 

We consider multi-matrix models where the dynamical variables are a set of M 
hcrmitian N x N matrices [Ai\%. Here i — 1, • • • ,M labels the matrices and 
a.b = 1, • • • , N are the 'color' row and column indices. The action S(A) for the 
matrix model is a U(N) invariant polynomial in S(A) — trS 1 Aj. 1 An example 
is the 2 matrix model 



S= tr 

Let 



-\[Ai,A 2 ] 



(1) 



dAe~ NS ^ (2) 

denote the partition function. The observables we are interested in are the 
correlation tensors Gi of the large N limit, the limit where N — > oo holding the 
coupling constants S 1 fixed. 



< > = < —A H ■ ■■ Ai >= - [ dAe- NS ^—A n ■ ■■ Ai 

1 Capital letters denote multi-indices S 1 = S' 11 "" 1 ", Aj = A il A i2 ■ ■ ■ A ip 
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Gh-i P = Jim < > (3) 

iv — >oo 

These are a complete set of observables in the N — > oo limit since expectation 
values of products of invariants factorize 

< $7!$^ • ■ • $/ r >=< >< $/ 2 > ■ • ■ < $/ r > (4) 

The Gi satisfy factorized loop equations (factorized Schwinger Dyson equa- 
tions): 

S Jllj2 G JlIj2 =5\^G I± G l2 (5) 

On the left side is an action dependent term while on the right is an anomalous 
universal term related to the non-commutative entropy explained in |3()j . 

It has been shown |17| that the integrals over matrices for the pure commu- 
tator squared action for a two matrix model are not convergent. To see this, 
consider the partition function, and go to the basis in which Ai is diagonal. 
In this basis, the integrand is independent of the diagonal elements of A2, and 
therefore diverges. The divergence is even worse if we consider the expectation 
of the trace of a polynomial involving A2, since then the integrand grows for 
large values of the diagonal elements of A^. It is thus necessary for us to reg- 
ularize the action. The simplest possibility is to add a quadratic term which 
ensures convergence of the matrix integrals: 

S = \{Al+Al)-£[A u A 2 f (6) 

This is the model we focus on. There is another reason to consider this two 
matrix model. An important property of Yang-Mills theory (in any dimension) 
is that its action leads to factorized loop equations whose action-dependent term 
is a derivation of the shuffle product of correlation tensors 2 . Therefore we would 
like to study a two matrix model whose action shares this property. It can be 
shown that the most general quartic 3 M-matrix model with this property is 

2 This remark will be explained in detail in a forthcoming paper on the algebraic structure 

of factorized loop equations 1331 . 

3 There are polynomial interactions of higher order with the derivation property. 
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tr 



S lll2 Ai 



1 i-^-ii , Ai 2 ]Ai 



' (Ai x i 



A,,. 



Ai 3 i 2 i 1 i 4 Ai 3 i 1 i 2 i 4 ) 



(7) 



where S 1 are arbitrary cyclically symmetric tensors. We can use a GLm(C) 
change of basis Ai = T?Aj to reduce the action to the canonical form where the 
covariance S lJ i— > |<5*-. Under such a change of basis 5(A) = S 8l "" ! ™A.j 1 ...j n i— > 



Thus the 



S(A) = S^-^T]l---T>A ix ... in and C, .... 
correlations Gi can be obtained from those of the canonical action. In the two 
matrix case (M = 2) the action reduces to 



S= tr 



±(A\ + Al)-»-[A x ,A 2] 



(8) 



There is only one independent coupling constant, the ratio of the coefficients of 
the quadratic and quartic terms. It is more convenient to study 



S 



tr 



-(Al + A 



[Ai,A 2 f 



(9) 



2 V 1 4, 4 L 

since it allows us to consider the pure commutator squared model in the m 2 — » 
limit. In the large AT limit, all correlations are functions of the single coupling 
constant m 2 . In the pure commutator squared model, the coupling constant 
is an overall factor in the action and can be scaled out, the dependence of 
correlations on it can be determined by dimensional analysis. By contrast, in 
the model we study, the coupling constant dependence of the correlations is 
to be dynamically determined, making it more analogous to Yang-Mills theory 
than the large N reduced models of M-theory [HI IT51 HI ITB1 fTBl fT7| 

Let us now summarize how we measure correlations by Monte-Carlo simu- 
lation. 



3 Monte Carlo Measurement of Correlations 

To measure the correlations numerically, we generate an ensemble of matrix 
configurations A^\ k = 1, ■ • • , n such that as n — > oo, matrix elements [Ai]% 
picked at random from this ensemble are distributed according to j^e~ NS ( A \ 
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The Metropolis algorithm |34l 1351 131)1 137| is used to create such an ensemble. 
We begin with a configuration A^p such as the zero matrices. The matrix 
elements are updated sequentially preserving hermiticity. At the k th time step, 
a candidate configuration is generated Bi = A\ +Wi + W}. Here for each i, 
Wi is a random N x N Weyl matrix. The only non-vanishing matrix element 
of Wi is picked at random from a uniform distribution of complex numbers in 
a square whose diagonally opposite vertices are at —A — V - 1A, A + \J— 1A. 

If B has a greater Boltzmann weight than A^ k \ the candidate is accepted, 

Af +1) = Bi if e -»WW-a{A™)) > 1 (10) 

and when B has a lesser Boltzmann weight than A^ k \ the change is accepted 
with probability e~ N ( SiyB ^~ s ( A{k) ^ and rejected otherwise 



Af +1) = Bi if 1 > e -^s(B)-s(A^)) > r 

Af +1) = 4 (fc) otherwise (11) 

Here r is a random number uniformly distributed in (0, 1). 

We can regard each matrix element as a continuous spin variable. One 
Metropolis sweep corresponds to \{N 2 — N) + N time steps of sequentially 
updating each independent matrix element of the M matrices. We perform 
a large number n s of such Metropolis sweeps, generating an ensemble with 
n = t;7i s (N 2 + N) configurations (n = 18000 for the measurements presented in 
this paper). 

It is shown |35l I36| that as n — > oo this algorithm produces a Boltzmann 
ensemble of configurations. The idea is that the above rules for making a tran- 
sition can be used to define a Markov matrix on the space of ensembles. Its 
eigenvalue of maximum modulus (=1) corresponds to the eigenvector labelling 
the Boltzmann ensemble. Thus one defines a contraction mapping on the space 
of ensembles whose unique fixed point is the Boltzmann ensemble. 

Once we have this ensemble, the correlations are computed: 

k=l 
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To extract the large N limit of the correlations, we do the measurement for 
several values of N = 10, 7, 15 and fit the results to the known N dependence 
for large N: 

• ^f-l- ' ■(•■■: ■ ^ (13) 

and extract the value of Gi 1 ...i p . 

We mention the main sources of error in these measurements. First, there is 
the statistical error (0( w "^? ce )) of truncating the ensemble of configurations at 
a finite value of n. This is estimated by the bootstrap [HHj procedure. Next there 
is the systematic error that could arise from the choice of initial configuration of 
matrices. This is estimated by changing the initial configuration slightly. The 
truncation (A) of the region in the complex plane from where the increments 
to random matrix elements are picked is a third source of error. As a practical 
matter we pick A so that the acceptance of the algorithm to candidate config- 
urations is roughly 50%. For sufficiently large n, we find the measurements to 
be insensitive to small changes in the value of A. The value of A needs to be 
increased in order to measure correlations of very high order accurately while 
holding other parameters fixed. Finally, there is the error in extracting the large 
N limit of the correlations from finite N data. 



4 Variational Principle for Large N Matrix Mod- 
els 

Let us summarize the variational principle introduced in jHOj . Given an action 
S(A) = tr S 1 Aj, for an M matrix model we want to determine the correlations 
Gj in the large TV limit. We found a variational principle fi(G) = x(G) — 
S I Gi whose extremization leads to the factorized loop equations S Jlij2 Gj 1 ij 2 = 
S I I llI ' 2 Gi 1 Gi 2 . The variation of S I Gi gives the action dependent term on the 
left while the variation of x gives the universal term on the right. The main 
difficulty was that the non-commutative entropy x is not a power series in the 
Gj. The space of correlations is a coset space of the automorphism group of 
the free algebra by the subgroup of measure preserving automorphisms. So we 
expressed x as an invariant power series on the larger space of automorphisms of 
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the free algebra, x is actually a non-trivial 1-cocycle of this group. The formula 
for x is given in |30| . we will only use a special case of it in this paper. 

Thus, to determine the correlations Gj, we must maximize the entropy x 
while holding the correlations Gi conjugate to the coupling tensors S 1 in the 
given action fixed. In other words, in the factorized loop equations, S 1 are the 
Lagrange multipliers enforcing these constraints. 

The resulting maximum value Xmax has a simple physical meaning. Let 
Sq(A) = i tr Y]fL x AiAi be the canonical gaussian action for an M matrix 
model. Suppose Z and Zq are the partition functions of S and Sq. Then 



(14) 



' Z ' 

i.e. —Xmax is the free energy or vacuum energy measured with respect to the 
canonical gaussian matrix model. 

Maximizing Q exactly is equivalent to solving the factorized loop equations 
exactly, which has not been possible in general. To determine the correlations 
approximately, we maximize Q on a subspace of the configuration space. The 
simplest subspace is that corresponding to the correlations of the multivariate 
Wigner distribution. This is the wignerian variational ansatz. On this subspace 

to[Gij] = ~ logdet[G y ] - S'd (15) 

where Gi are the correlations of the multi-variate Wigner distribution. They are 
determined in terms of the two point correlation matrix by the planar analogue 
of Wick's theorem: Gijki — GijGki + GuGjk, etc. Thus we regard the matrix 
elements of Gij as variational parameters and maximize Q for the given action S. 
To summarize, the wignerian variational ansatz gives the wignerian correlations 
that best approximate the true correlations of a given matrix model. Best 
approximation here means the one that maximizes entropy while holding the 
correlations conjugate to S 1 fixed. 

We now compare our wignerian variational ansatz, first with the exact so- 
lution of a two matrix model studied by Mehta and then with the numerical 
solution of the gaussian + Yang- Mills two matrix model. 
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5 Two Matrix Model studied by Mehta 

We recall here (for details see [30]) the comparison of our variational ansatz 
with Mehta's exact solution of a two matrix model. In Mehta finds the 
exact vacuum energy of the two matrix model with action 



S(A, B) = tr i(A 2 + B 2 - cAB - cBA) + |(,4 4 + B 4 ) 
From his solution, we may extract the exact values of 



E .. (9 , c) = _ lta _L log |||, Gi - BmdGlw ( i„ 

in the strong and weak coupling regimes. These are compared with our varia- 
tional estimates below. For small g and c = \: 



E ex (g, 
E var (g, 



GTb{9,\ 



Gaaaa(9i 

^AAAAKSi 



AAA +1.780- 8.7V 
.144 + 3.56g- 23.7g 2 



-- 4.74.g + 53.33 5 2 

- - 4.74.g + 48.46s 2 
3 



32 

y 

32 



34.960 H 

31.610 + 368.O20 2 



(18) 



For strong coupling and arbitrary c: 



E ex (g,c 
E var (g,c 

G e A x B (9,c 
G v A a B(9,c 



Gaaaa(9i c 



G v aaaa^9i c 



1 1 „ 3 
-log0+ -log3- - +••■ 

ilo g + ilo g 2 + ^ + O(i) 

2 2 V 8 3 9 



as — * oo 
" " " ^ 



20 (2g)i 



9" 



(20)1 



0{- 



(19) 
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We see that both for strong and weak coupling, our wignerian variational 
ansatz provides good estimates for the partition function and correlations. We 
now consider the gaussian + Yang Mills two matrix model, where we compare 
our variational estimates with Monte-Cairo measurements. 



6 Gaussian + Yang-Mills two matrix model 

Let us specialize to the two matrix model with action (m 2 > 0) 



S(A) = tr 
For the wignerian ansatz, 



(20) 



1 m 2 1 

Q[G] = - logdet[G y ] - — (G u + G 22 ) + -(G 1212 - G 1221 ) (21) 

Due to the A\ «-> A 2 symmetry of the action we can assume a variational matrix 



Gi 



a (3 
P a 



Since < j?A\ > > and < jj(Ai - A 2 ) 2 > > 0, we must maximize 



Q(a, 0) = \ log(a 2 - 1 ) ~ m 2 a + i(/3 2 - a 2 ) 



in the region a > and a > (3. We get 



(22) 



(23) 



G ll = G 22 =a=\ [ jl +■—-■—, G 12 = G 21 =j3 = (24) 

Figures^and|21compare the variational two point correlations with Monte-Carlo 
measurements for a range of values of m 2 . 

All other correlations can be expressed in terms of these. For example, the 4- 
point correlations are (the rest are determined by cyclic symmetry and A\ <-> A 2 
exchange symmetry) 



Gun — 2a z \ G\ 2 \ 2 — 0; Gi 22 i — a 2 ; Gm 2 — 



(25) 
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Figure 1: log 10 [Gn] versus log 10 [to 2 ]. Solid line is variational estimate, dots 
are the Monte-Carlo measurements. The approximation becomes poor for small 
values of to 2 
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Figure 2: G12 versus log 10 [m 2 ] 



Log Gi 




2 J Lo 9 m 



Figure 3: log 10 [Gnu] versus log 10 [m 2 ] 
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Figure 4: G1212 versus log 10 [to 2 ] 
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Figure 5: log 10 [Gix 22 \ versus log 10 [m 2 ] 

Figures |3 0] and [5] compare variational estimates (solid lines) and Monte-Carlo 
measurements (dots) of Gun, G1212 and G1122 for 10~ 3 < m 2 < 10 3 . The n 
point pure A\ (or A2) correlation is given by the Catalan numbers 



Gm...i — G222-2 
More generally, 



c%a 2 = fsAUWTT^\ a 2 , if n is even; 



G M = { -ItTKlW-' — (26) 



0, if n is odd. 



Gn...i22---2 = G( ril )(„ 2 ) — G( ni )G(„ 2 ) 



Gi 



l-"122"-2H".122"-2 



G( T i 1 )(„ 2 )( Il3 )(„ 4 ) — G(„ 1+ „ 3 )G(„ 2 )G(„ 4 ) 
+G( ni )G(„ 3 )G(„ 2+ri4 ) — G( ni )G(„ 2 )G(„ 3 )G(„ 4 )(27) 



We mention these since they are useful in estimating expectation values of Wil- 
son loop- like operators. The variational estimate for vacuum energy is 



E V ar(g) = - Jim ^log 



z(g) 
Z(0) 



= — log a = — log 



m 
T 



(28) 



6.1 Wilson Loop Operators 



It is also interesting to see what the wignerian ansatz says about the 2-matrix 
analogue of the expectation of the Wilson loop in the large TV limit. 
Wilson Line: The simplest analogue is a 'Wilson line' the analogue of the 
parallel transport along a line of length I in the A\ direction 



tr 



W Une {l)= lim < ^e lMl > 



N 



(29) 
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Figure 6: Wu ne (l) for m 2 = 1. Dots are numerical and solid line variational 
estimate. 

For the wignerian ansatz, we get (using eqn (|26Jl . J n (z) is the Bessel function 
of the first kind) 



oo / -j\2k 

71=0 '' 



-.Ji(2ly/a) ~ 3- cos — 2l\fa) as I — > 00 (30) 



Wu ne (l) is a real- valued function of real I since the odd order correlations vanish. 

Thus, for the wignerian ansatz, the expectation value of the 'Wilson line' is 

oscillatory but decays as a power Z -3 / 2 . For small I, Wn ne (l) — > 1 — \ai 2 + 
2 ,4 _ 

■ • •. Figure [6] compares this ansatz with Monte-Carlo measurements for 

to 2 = 1. The behavior both for small and large values of I is well captured by 
our ansatz. 

L shaped Wilson Line: For an L shaped curve, we define 

W L (l) = lim < ^- e ltAl e aA2 > (31) 

For the wignerian ansatz (use eq. l|27)l: iF 2 (a, h;z) — YlnLo (hi)°(6 2 ) nT * s a 
generalized Hypergeometric function with (a)„ the Pochhammcr symbol), 



W L (l) 



E 

ni ,712 =0 
00 

E 

fci,fc 2 =0 



W , , lim < ^U?^ > 



(-lo) 



2„,\fci+fc 2 



fci!(fci + l)!fc 2 !(fe 2 + 1)! 
= Vf-i 2 af - 1 2> 

^ 0Fr(n + i)r(7i + 2)r( ? i + 3) 
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Figure 7: Wl{1) for m 2 = 1 



W L (l) 



1 F 2 (-;{2,3};-4/ 2 a) 



(32) 



5a' 



As before Wl(0 is real for real L For small i, Wl(Z) — > 1 — al 2 
This is compared with the numerical calculation in fig. 0for m 2 = 1. Both the 
small / behavior and decay for large I are captured by our estimate. 
Wilson Square: The analogue of the parallel transport around a square of 
side I in the A\ — A2 plane is 



tr 



W square (l)= lim < ^ e ilA le UA 2e -ilA le -ilA 2 



> 



(33) 



N-^oo N 

In the wignerian variational approximation, W squar e{l) is real- valued since odd 
order correlations vanish. Using ea. (|27|l we get 



W square (l) = ^H 2 a)"T!(n) + 2(l 2 a) ^(-/ 2 a )«T 2 (r 



(34) 



n=0 



n=0 



where, 



Tx{n) = Y, 



2c kl+k3 Ck 2 c ki - n^i£fcj 
ILf_J2ki)\ 

fci>0,feiH \-k 4 =n ' 

Ck 1 +k 3 + lCk 2 Ck 1 



rp I \ ST^ L fci+fc 3 + l^fc 2 L fc4 /ar\ 

^ fel ,o^. +fc4 ^(^ + 1 ) ! ( 2fc2 ) ! ( 2fc 3 + l)!(2fc 4 )! 1 j 



For small i, W sguare (i) — * 1 — I 4 a 2 + • • •. The expectation value of the 

Wilson loop is a rapidly decaying function for large values of I. It would be 
interesting to find the asymptotic rate of decay. It is oscillatory but a posi- 
tive function, unlike Wu ne . These variational predictions are confirmed by the 
numerics, in fig. [H] for m 2 — 1. 
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(l) for m 2 = 1 



The variational ansatz does a very good job of estimating the Wilson loop 
averages over the entire range of values of I studied, for m 2 of order unity or 
more. 

7 Summary and Discussion 

We find that despite its simplicity, the wignerian ansatz for our entropic varia- 
tional principle is remarkably good at estimating correlations. For the exactly 
solved model studied by Mehta, our ansatz works well both for strong and weak 
coupling. For the gaussian + Yang-Mills two matrix model, our estimates for 
correlations and expectation values of Wilson loop operators are accurate for 
moderate and weak coupling when compared with Monte-Carlo measurements. 
When the commutator squared term dominates the action, the wignerian ansatz 
becomes poorer as an approximation. This is to be expected. The measured 
correlations grow without bound as we approach the pure commutator squared 
limit. This reflects the divergence of the matrix integrals for the pure commu- 
tator squared interaction. 

We can go beyond the wignerian ansatz for the entropic variational principle. 
In particular, this will allow us to improve on the wignerian ansatz and also 
estimate correlations that are identically zero for the wignerian ansatz. We will 
address this question in a future paper. 

In another direction, we hope to extend our approximation methods to su- 
persymmetric matrix models in order to make predictions about the matrix 
models of M-theory and superstring theory mentioned in the introduction. 
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A conceptual shortcoming of the numerical procedure used in this paper is 
that we measured correlations for several finite values of N before extrapolating 
to the N — oo limit. Is there some way of computationally determining the 
N = oo correlations directly? This is a challenging problem that would likely 
require new ideas from non-commutative algebra, geometry, probability theory 
and computer science to determine directly the N = oo correlation tensors 
without integrating over matrix elements. 

It would be interesting to make precise the field theoretic connection between 
general multi-matrix models and supersymmetric gauge theories with adjoint 
chiral superfields. Our methods for estimating the correlations of the former 
can potentially shed light on the expectation values of operators in the chiral 
ring of the supersymmetric gauge theory. 

Our investigations have focussed on matrix models with a finite number of 
matrices. It would be desirable to have a similar theoretical approach based 
on approximation methods over and above the large N limit, for matrix field 
theories (beyond c = 1 matrix quantum mechanics) such as Yang-Mills theory. 
This would complement the lattice gauge theory Monte-Carlo efforts (see eg. 

my 
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